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ABSTRACT 

The selfconsistent steady state solution for a strong shock, significantly 
modified by accelerated particles is obtained on the level of a kinetic description, 
assuming Bohm-type diffusion. The original problem that is commonly 
formulated in terms of the diffusion-convection equation for the distribution 
function of energetic particles, coupled with the thermal plasma through the 
momentum ffux continuity equation, is reduced to a nonlinear integral equation 
in one variable. The solution of this equation provides selfconsistently both the 
particle spectrum and the structure of the hydrodynamic ffow. A critical system 
parameter governing the acceleration process is found to be A = M^^/'^Ai, 
where Ai = rjpi/mc, with a suitably normalized injection rate r), the Mach 
number M ^ 1, and the cut-off momentum pi. 

We are able to confirm in principle the often quoted hydrodynamic prediction 
of three different solutions. We particularly focus on the most efficient of 
these solutions, in which almost all the energy of the fiow is converted into 
a few energetic particles. It was found that (i) for this efficient solution (or, 
equivalently, for multiple solutions) to exist, the parameter C = r)y/popi/mc must 
exceed a critical value (ci ~ 1 {po is some point in momentum space separating 
accelerated particles from the thermal plasma), and M must also be rather large 
(ii) somewhat surprisingly, there is also an upper limit to this parameter (iii) the 
total shock compression ratio r increases with M and saturates at a level that 
scales as r oc Ai (iv) despite the fact that r can markedly exceed r = 7 (as for 
a purely thermal ultra- relativistic gas), the downstream power-law spectrum 
turns out to have the universal index g = 3 Y 2 over a broad momentum range. 
This coincides formally with the test particle result for a shock of r = 7 (v) 
completely smooth shock transitions do not appear in the steady state kinetic 
description. A finite subshock always remains. It is even very strong, rg ~ 4 for 
A <^ 1, and it can be reduced noticeably if A ;^ 1. 

Subject headings: acceleration of particles, cosmic rays, diffusion, hydrodynamics, 
shock waves, supernova remnants 
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1. Introduction 



Strong astrophysical shocks are widely believed to be the sites where the cosmic rays 
(CRs) are born. Although the test particle calculation of the CR spectrum is straightforward 
(Krimsky 1977; Axford, Leer & Skadron 1977; Bell 1978; Blandford & Ostriker 1978), 
the acceleration efficiency, in other words the fraction of incoming flow energy that is 
converted into CR-gas internal energy, cannot be obtained within this theory. There are 
two major problems associated with diffusive shock acceleration. First of all even in the 
simplest test particle theory a self-similar power-law solution, isotropic to lowest order, 
is possible only at sufficiently large momenta where the solution is independent of any 
momentum scale and depends only on the shock compression ratio. Thus, some difficulty 
in this description occurs already at lower energies where the thermal and the bulk flow 
velocity may enter the solution. Consequently, the amplitude of this power-law spectrum is 
virtually unknown. It is therefore necessary to calculate, flrst of all, the injection rate, i.e. 
the number of particles that feed the acceleration process. Injection is believed to operate 
at the plasma subshock. Surprisingly, until recently there were no attempts to attack the 
injection problem analytically. At the same time this problem was studied numerically 
( lEUison 19851 ; [Quest 1988|; [Scholer, Trattner fc Kucharek 1992|; |Kang fc Jones 1995|) and 



phenomenologically ([Lee 1982 ; Zank, Webb fc Donohue 1993 ) in great detail and from 



different points of view. Given the subshock strength and in terms of the parameterized 
distribution of thermal particles leaking into the upstream region from downstream, the 
injection rate has been calculated analytically by |Malkov fc Volk (1995), (MV95 thereafter' 



The distribution of these leaking (and/or directly from the shock front reflected) ions that is 
crucial for the injection theory is in turn a problem of its own in coUisionless shock physics 
( ^agdeev 1966| , [Kennel fc Sagdeev 1967| ). It almost certainly depends on shock parameters 
like the orientation of the magnetic field, the ratio of upstream thermal and magnetic 
pressures, the Mach number, etc. One simplified formulation of this formidable task that 
also allows one to calculate this distribution in a closed form, suitable for injection, has 
been suggested by this author ( [Malkov 1996| , M96 hereafter) for a parallel shock, where the 



direction of the shock normal and the magnetic field coincide. For other types of shocks 
the question is unsolved. It may be said then that the injection problem in general has 
become more a problem of the shock dissipation mechanism, rather than the problem of 
shock acceleration. It is however also true as we shall see in the sequel that these two latter 
aspects of coUisionless shock theory cannot be treated independently. 

The second major problem, the impact of the accelerated high-energy particles on the 
acceleration process itself may be even more dramatic. While the injection rate can be 
calculated given the subshock conditions, even though subject to the above hmitations, 
and the actual subshock strength can be determined afterwards (when the pressure of 
accelerated particles is finally obtained, see Appendix A), the calculation of the spectrum of 
dynamically important high-energy particles should be intrinsically nonlinear. The reason 
is that the energetic particles can significantly modify the fiow structure over a large spatial 
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scale ~ K,{pi)/ui, where K,{p) is the momentum dependent diffusion coefficient (in the Bohm 
hmit one has k oc p). Here pi is the upper cut-off momentum and Ui is the ffow speed 
far upstream. Hence particles with different energies "see" different shock compression. 
Moreover, they produce the gradient of the flow velocity by themselves, and therefore couple 
the length scale with the momentum scale. This means that a momentum scale-free particle 
spectrum is no longer possible unless the velocity profile is also scale invariant. Since 
in a strongly modified high Mach number shock the total compression ratio exceeds the 
conventional value of 'four' markedly, the partial pressure of stationary accelerated particles 
becomes a nonintegrable function of momentum without an upper momentum cut-off. 
Therefore, in reality the steady state acceleration along with the underlying flow structure 
will critically depend on the losses at the upper cut-off momentum pi or, more precisely, 
they will be determined by the balance between these losses and the injection around some 
slightly suprathermal [| momenta p ^ Po^ Pth (see appendix A). Moreover, these losses have 
the effect of increasing the total compression ratio, boosting CR production even further. 
One sees that also the form of the spectrum obtained within the test particle theory is 
dubious since the backreaction of accelerated particles on the flow is significant. 

In this paper we present a solution of the following problem. Consider a strong, 
stationary, plane shock, propagating at the Mach number M ^ 1. Suprathermal particles 
are steadily injected at the subshock and are then partly accelerated up to the cut-off 
momentum pi ^ mc. To be determined (as functions of M and pi) are: (i) The spectrum 
of accelerated particles, (ii) The flow profile across the shock u{x). The latter implies the 
total compression ratio and thus the acceleration efficiency, again as function of M and pi. 
Of particular interest are the question of the uniqueness of the solution, and the strength of 
the subshock. 

Physically, the injection rate t] does not belong to the input parameters which are 
merely M and pi. It can be calculated selfconsistently with the help of (MV95, M96, 
and |Malkov fc Volk 1997| , (MV97)), given the subshock strength to be obtained below. 
Nevertheless, under certain restrictions explained in Appendix A, injection can be treated 
independently of the large-scale shock modification studied in this paper. We will pursue 
this approach in what follows, considering r] as another input parameter. 

There is no necessity to stress that this paper is not the first one to attack the 
problem of nonlinear shock acceleration analytically. Besides the well known two-fluid and 
three-fluid approximations introduced by Axford, Leer, & Skadron (1977), Drury & Volk 
(1981), and, including the scattering wave field, by [McKenzie fc Volk 198^ on the one hand, 
and the perturbative kinetic studies performed by Blandford (1980) and Heavens (1983) on 
the other, there exist a few separable solutions based on rather special assumptions about 



^We regard particles downstream as suprathermal if they can be scattered upstream 
elastically. 
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the functional dependence of k{p,x) {e.g., Drury, Axford & Summers 1982; |Webb et~al. 
1985|) . An alternative treatment was suggested by Eichler (1979, 1984). The key step of 
his approach consists in replacing the true solution of the diffusion-convection equation ad 
hoc by a Heaviside function which is coordinate independent up to some distance xo{p) 
upstream and is zero beyond it. Unfortunately, such a 'weak solution' does not satisfy 
the convection diffusion equation which may be proven by substitution. Physically there 
is indeed no other scale height upstream for the energetic particles than the diffusion 
length k{p)/u which could justify the introduction of a step-like behavior of the particle 
distribution. The width of the transition zone between the spectrum downstream and its 
far upstream (zero) limit increases with p exactly as does the length parameter xq {p) where 
the spectrum supposedly vanishes abruptly in Eichler's model. As it was pointed out by 
Eichler (1984), such a behaviour of the solution would be possible if k{p,x) were extremely 
large for x < Xq{p) and zero otherwise. Therefore, this solution belongs also to the category 
of solutions that require very special form of diffusion coefficient k. There is no physical 
reason, however, to ascribe such a property to the CR diffusivity. 

Notwithstanding the above criticism, the approach developed by Eichler contains the 
fruitful idea of 'algebraization' of the problem. In other words, it offers a way to reduce 
the number of independent variables by the introduction of some functional link between 
the X and p variables through a special boundary x = xo{p) in phase space. In addition, 
the issue of interrelation between injection and losses and their critical influence on the 
shock formation and acceleration efficiency was first addressed by Eichler (1979, 1984). 
In principle some form of algebraization seems to be necessary for a successful analytic 
treatment. Otherwise computer studies will remain the only possibility to explore this 
difficult problem. In contrast to Eichler we derive an algebraization that is based on the 
powerful technique of integral transforms. 

We shall focus predominantly on the so called efficient solution (Fig.l), in which almost 
all the shock energy is converted into CR-gas internal energy. Its relation to the inefficient, 
test-particle solution is the subject of a companion paper (Paper II). By an extension of 
the technique developed below, we present in Paper II a unified description of all three 
solutions along with the corresponding bifurcation analysis. 

The outline of this paper is as follows. In sec. 2 we introduce the necessary equations 
and discuss some of the basic approximations used in the further analysis. In sec. 3 we find 
a formal solution for the particle distribution in the shock precursor, that depends on the 
unknown postshock spectrum and on the hydrodynamic velocity in the shock transition. We 
introduce a spectral function that couples the thermal gas with the energetic particles and 
plays a central role in the analysis. In sec. 4 we derive equations for the spectral function 
and for the hydrodynamic flow, investigate analytic properties of the spectral function and 
solve the equation for it. Sec. 5 deals with the flow profile. In Sec. 6 the particle spectrum is 
restored. We summarize and discuss the results in sec. 7. 
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2. Basic equations and assumptions 

Consider a strong CR-modified shock propagating in tlie positive x-direction in a cold 
medium with the gas kinetic pressure Pgi and the mass density pi. In the shock frame of 
reference the steady mass flow profile is defined as follows 



U{x) = l (1] 




where U2 is the (constant) downstream mass velocity, ^(0+) = uq, and u{x) — > mi = const, 
as X — oo. The isotropic part of the particle distribution at sufficiently high momenta, and 
for all X, is governed by the well known diffusion-convection equation (Parker 1965; Gleeson 
& Axford 1967) 

— (ug - K(p)—\ = -—p^ (2) 
dx \ dx J 3 dx dp 

Here the number density of CRs is normalized to Angdp/p, k denotes the particle diffusion 
coefficient that is assumed to be a monotonically increasing function of momentum p. To 
be specific, we will assume that K,{p) = i^op/po = f^ip/pi- For x < the solution to eq.(0) is 
taken to be U{x) = —U2 = const and g = go{p)- For x > g{x,p) satisfies the equation 

d f dg\ Idu dg 

dx \^ ^ dx) Sdx^dp 

which should be solved together with the conditions of overall mass and momentum 
conservation 



pu = piui, (4) 

Pc + Pini+Pu' + P^li-Y = Piul + P^l (5) 



u 



Here p(x) is the mass density, pi = p(oo), 7 is the specific heat ratio of the thermal plasma, 
Pc is the CR pressure 

47r ^ 2 [^^ Pdp 



Pc{x) = —mc / —==g{p,x) (6) 
3 Jpo \Jp^ + 1 

and no seed particles are present, i.e. Pc(oo) = 0. The particle momentum p is normalized 
to mc. We regard as CRs all the particles that occupy the region in momentum space 
between the injection momentum p^ and the cut-off momentum pi, and we assume that 
Pq < \ P\. In Bernoulli's integral (^) we have isolated the contribution of suprathermal 
particles Pinj that do not belong to the thermal upstream gas and should rather be regarded 
as a low-energy part of the CRs. These particles appear just upstream from the gaseous 
subshock as a product of thermal leakage from the downstream medium and have practically 
nothing to do with the adiabatically compressed upstream thermal plasma represented 
by the last term on the l.h.s. in eq.(^ (see Appendix A). In view of the astrophysical 
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significance of strong shocks, we focus on the high Mach number hmit. For simphcity we 
confine our consideration to the case in which not only M ^ 1 but also 

where = piuf/jPgi. Thus eq.(|^) rewrites 

U{X) + ^ (Pe + Pinj) = Ui (8) 
PlUi 

One remark should be made concerning the approximation in eq.(^. Namely, we neglected 
the gas pressure by virtue of assumption (0), comparing it with the ram pressure piuf. 
If the shock is significantly modified by CRs, one should compare the gas pressure with 
piuf — Pc instead and the condition may become more stringent, uui/uq <C 1. However, 
this may potentially be important only close to the subshock where also the Pinj-term may 
be important, and the problem can hardly be resolved without turning to the subshock 
internal structure itself. We will specify the corresponding length scale and discuss the 
validity of eq.(||) in Appendix C, after the flow structure is determined. The last equation 
we need is that for the gaseous subshock, in which equation we neglect the energy losses 
from the gas due to injection: 

^ = ^''^ , (9) 

It is important to note that the subshock can still be diminished significantly due to the 
CR pressure and at least formally can even be smeared out completely when u = 2mo/mi, 
without violating our assumption (0). To summarize, eqs.(||,0,|), and (H) form the basis for 
the further analysis. 



3. Approximate solution to the diffusion-convection equation 

The main purpose of this section is to gain enough information about the solution of 
the diffusion-convection equation without specifying the flow profile u{x). To this end 
we will use the flow potential 

^ = r udx (10) 



as a new independent spatial variable instead of x. But first, we transform eq.(H) to the 
following integro-differential equation 



gix,p) = e-*(^')/'^h(p)- 

f) roo rill 

(11) 



^ r e^^^'y^dx'p^ rp^g{x",p)dx" 
6k Jo op Jx' ax ' 
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K — 

dx 



Once the solution of this equation is found in terms of the unknown downstream spectrum 
go{p), the latter can also be determined from the equation immediately following from 
eq.(l) 

One sees that iterates g^""^ of eq.(|ll]), starting with 

g^''^ =go{p)exp{-^{x)/K) (13) 

yield a uniformly valid expansion in the case {ui — uq)/ui ^ 1 (weak modification). In 
fact, is already a good approximation in this case. Let us start our consideration of 
the strongly modified shock from the region <C 1. Substituting g^^^ into the r.h.s. of 
eq.(|llD we can generate a formal Neuman series in Indeed, calculating the internal 

integral by parts and retaining formally only the leading term in !/«;, as a first iteration we 
get 

= g,e-^'^ (l - (14) 

where 

3 amp 

To the same order in we can rewrite g'^^^ in the same form as g^^\ but with a 
renormalized diffusion coefficient 

^«=^o(p)exp|-i±^M/| (16) 

Since (5 may be rather small numerically, {e.g., for the test particle solution with the 
compression ratio r = 7, one obtains [3 = 1/6) and because g^^^ provides a correct 
asymptotic behavior of the solution of eq.(|^) for \I' — > oo {du/dx — > oo), one may conjecture 
that the solution in the form of eq.([T6|) is also good for larger values of \E'/k. That this is 
true, is demonstrated in Appendix B. Therefore, we may take 

<7 = <7o(p)exp|-i±^vI/| (17) 

as a formal solution to the diffusion-convection equation. The function go{p), or 
equivalently, P is yet to be obtained from eq.(|r2p whereas the function \l/(x) must come 
from the Bernoulli's integral (Mj. 



3.1. Equation for the particle spectrum 

The solution (|l^) is valid for < 6^^, where 5 <^ 1, and its region of validity has, 

in fact, also a lower bound given by eq.(|B3|) (see Appendix B). As it will be shown later. 
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this restriction concerns only a very small region at the origin for the solutions of interest. 
Clearly, this may be important for dg/d'^ at small \1/ but not for g itself since the latter 
is continuous at \1/ = 0. Thus, some care should be exercised in calculating the l.h.s. of 
eq.(|l2D. In fact, it is sufficient to rewrite this equation with the help of eq.(p|) as 



'"v-^^r?^?(-,pM- (18) 



d\np Au Augo d\np Jo+ dx 

This equation can also be derived directly from eq.(|^) by means of integration over x 
between 0— and +oo (Eichler 1979). In contrast to eq.(|T2p it does not contain the uncertain 
spatial derivative of the solution at the origin and we may substitute the solution (|1^) into 
eq.(|l8D which yields 

1 dV 

(19) 



3 d\np Au + V 

Here we have introduced the function V that will play a central role in our further analysis 



V(p) = [ ^dxexp 
Jo+ dx 



--^(x) 



K 



(20) 



This function reflects explicitly a degree of shock modification. In an unmodified shock 
V{p) = 0, since then du/dx = in the upstream region; the spectral index g = 3/5 is just 
the conventional q = 3m2/(mi — M2), (see eq.(|l^)). In general < V{p) < ui — uq and 
V{p) — > Ml — as p — > 00. Even if the shock is appreciably modified, one may show 
(Appendix C) that at small p Po 

V{p) < Am (21) 

and V may be neglected in eq.(|TP|). The spectral index then corresponds simply to the 
subshock compression ratio, and at lower momenta we have 

3u, / - ^ 



qc::iq, = go{p) = Qinj - (22) 

Uo - U2 \Po J 

The injection solution (Appendix A) produces essentially the same asymptotic result for 
p ^ po, yielding thus the injection rate Qinj. The solution go{p) can be obtained then for all 
p using eq.(^). To this end an independent equation for V{p) should be derived. This will 
be done in the next section. 



4. Integral equation for the spectral function 

In the previous section we have seen that the function V{p), that we term the spectral 
function, is directly related to the spectral slope through eq.(|l^). If V{p) were defined as 
an integral starting at 0—, the subshock jump would also be incorporated into the spectral 
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function as a point set component of the measure u{x) in the integral (^). We will use this 
'full' spectral function occasionally, denoting it by V{p) = V{p) + Am. We emphasize that 
in addition to the link with the particle spectrum, this function contains all the information 
about the hydrodynamic flow structure. Indeed, introducing a new variable 



k{p) 



(23) 



the function V can be recast as the Laplace-Stieltjes transform of the flow velocity m(\1/) 



V{s) 



e-'^'dui^), and V{s) 



(24) 



that maps 'u(\l') onto V{s), or in other words, V provides the spectral representation of 
the flow m(\I/) in terms of the particle momentum p by virtue of eqs. (|23|j2^ ). Yet another 



interpretation of V{p) is that a particle with momentum p samples the flow compression 
6u = u{x) — uq upstream, where x is defined by the relation (1 + /5)\1/(x)/k(p) ~ 1. 

After these short remarks concerning the usefulness of the spectral function V for 
describing the CR-gas coupling, we note that the only equation available to derive an 
equation for V is the Bernoulli's integral (|^). Using eq.(|T^) the latter can be rewritten as 



— + F{^) 

ax 



Ui 



where 



F(^') = r]Ui 



PI pgodp 



po 



exp 



Here the injection rate is given by 



V 



Arc mc^Qin] 



pul 



k{p) 



go mricC 
3 piuj 



Pi 



(25) 



(26) 



(27) 



and ric ^ pi/m denotes the number density of the CRs. The low energy asymptotics of go 
is now fixed by = (po/p)^^ p Po, i.e go{p) = Q\n]9o{p) (see eq.(||)). We will omit the 
tilde sign in what follows. The quantity r] is regarded as a known function of the subshock 
parameters, that can be inferred from the injection theory MV95, M96, MV97. 

For \E' > the pressure integral (^) is cut by the exponent at small k{p) <^ \E'. This 
means that starting from some \1/ ^ k(po) the first term on the r.h.s. of eq.(^) becomes 
virtually independent of po since, due to the normalization (?o(po) = 1, the function go 
depends also on po in such a way that rjgo is po-invariant. For smaller \1/ (0 < \1/ < k(po)) 
only the sum of the two terms on the left hand side (l.h.s.) is independent of po and the 
contribution Pinj that has the same scale in \l/ might be important. This is, however, a very 
small region compared to the total precursor height \E'p ~ i^{pi) ^ i^{po), in the case of 
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strongly modified shocks we are interested in. Therefore we may ignore here this region and 



omit the term P\ 



my 



Taking the x derivative of eg. , multiplying the result by exp(— s\E') and integrating 
with respect to x we get 

roo 

V{s) + / F'(^)e-^*d^ = (28) 
Jo 

Returning to the momentum p (eq.(E3])), for efficient solutions from (E6|) we obtain 



V{p) — rjUi 



Pi p'dp'go{p') s{p') 
Vp''^ + 1 s{p') + s{p) 







(29) 



(Note that for inefficient, weakly modified solutions, the function s{p) should be set to 
s = 1/k in accordance with eq.(|l3D, instead of eq.(|^).) Resolving eq.(|l^) for and using 
the normalization go{po) = 1 in place of the last equation we obtain the following equation 
for V 



V{p) = TjUi 



Pi p'dp' s{p') Au + V{po) 



po 



^/p^TT S{p') + S{p) Au + V{p') 



exp 



d In p" 



Po 



Au + V{p") 



(30) 



where the function s may be written as (eqs.(|T5|,|T9| 



s{p) 



uo + V{p) + 



1 dV 
3 91np 



k{p) {Au + V{p)) 



(31) 



in the case of efficient solutions, and s = 1/k{p) for inefficient ones. Equation ([30| ) is a 
nonlinear integro-differential equation for the function V{p). Mathematically, the advantage 
of this equation is that it is only an equation for one function of one variable. By contrast, 
the original equations (^,^) form a quasi-linear system of the partial differential equation (0) 
in two variables x,p under the algebraic condition (§). Moreover, provided that eq.(|30D is 
solved, 



eqs. (^ , p8D and (|l^) allow us to calculate both the flow proflle u{x) and the particle 
distribution g{x,p). Though eq. (|30|) is fairly complex, it can be easily solved numerically. 
It might be tempting to do so. However, such a shortcut would not allow much further 
insight into the physical nature of the solutions. Therefore, in this paper, we shall pursue 
an analytic approach with the aid of a number of simplifying approximations. 



4.1. Simplification of the integral equation for the spectral function 

We flrst introduce a normalized spectral function 

_ A^H-VM , m (32) 

Au + Vo Au* ^ ' 
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where Vq = V{po), Au* = Au + Vq, and substitute the relativistic Bohm diffusion coefficient 
Using a new variable t = Kos{p) in place of p, eq.(]5^) rewrites 



. /■*" dt' 1 
Here we have introduced the notations 



3U2 f^o dt" 



OAu* Jt' t"J{t") 



+ 



Au 
Au* 



A = r]po 



Au*'' 







An' 



Ul = U0 + V0 + 



1 dV 



3 d\np 



P=PO 



Pi 



V{pi)+uo + - 



3 d\np 



p=pi. 



[V{pi) + Au 



-1 



^^^«1 
Pi 



(33) 



(34) 



(35) 



Several comments should be made on eq.(|33D. Through an argument similar to that 
followed by eq. (pH]) we have simplified in eq.(|33|) the contribution of the nonrelativistic 
particles {po <^ p < 1) and used the Bohm diffusion coefficient n oc p for all p instead of a 
more general form of it, e.g., k cx p^/ \/l + p'^. We thus choose the simplest realistic form 
of the diffusion coefficient. Clearly such a choice suggests logically to replace a/1 + p'^ — >• p 
in eq. (pUD . The most serious potential problem that might be caused by such a simplified 
treatment is to overlook a relativistic peak at p ~ 1. The situation here is quite similar 
to that considered in ( Malkov fc Volk 1996[) , where the contribution of this peak to the 
CR pressure is shown to be negligible in comparison with that of the region 1 <^ p < pi 
for the efficient solutions considered further in this paper. These are mainly technical 
arguments for such a simplified treatment of eq. (pOD . The physical argument, as before, is 
an insignificant dynamical role of nonrelativistic particles. Besides, we replaced ds/dp by 

ds , . 
— ^ — s[p) 
dp p 

which is strictly valid only for V{j)) ^ Uq. Except for the region p ^ po, this is precisely the 
property of efficient solutions, as we shall see. 



Our next simplification of eq.(P3D is based on some assumption about its solution and 
this assumption will be justified later. Since for efficient solutions J(ti) ^ ^(^o) = I7 J(t) 
is large in the dynamically most important part of the interval {ti,to), i.e. for t ^ to or, 
equivalently, p ^ po- If we assume that the estimate J > Ct^'^ holds, with some positive 
constants C and A, then t' in the exponent of eq.(^3l) can be replaced by t' = ti, since the 
integral will be dominated by the upper limit to. Thus eq.(p3) can be written as 



J(t) = A* 



*o dt' 



+ 



Au 



H t' + tt'J{t') Au 



(36) 
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where we introduced the constant A*, to be determined from the following solvability 
condition 

A* A /■*" dt 

A = A exp 



OAu* tJ{t) 



(37) 



It will be specified after the solution to eq.(|36|) is found. Now, it is convenient to rescale 
variables in eq.(|36D as follows 



t 



T 




J 



(38) 



where 

and then eq.(^) rewrites as 

F(r) 



1/^ dr' 



Au etn 



e t' + T t'F{t') Au* V A* 



(39) 



The goal of the next two subsections is to find an asymptotic (e <^ 1) solution of this 
equation. 



4.2. Analytic properties of the spectral function 

The kernel of eq.(p9|) belongs formally to the Carleman type (see, e.g., [Muskhelishvili 



1953|) and the theory of linear equations of this type is very well developed. The problem 
is that eq.(^) is nonlinear and it is also well known that the solutions in this case can be 
fundamentally different from those in the linear case. In particular, multiple solutions and 
bifurcations may occur. A simple iterative approach would be quite reasonable for seeking 
inefficient solutions of eq.(PDD, or better, of the more general equation (pUD. Q In contrast 
to inefficient solutions, the efficient solutions being strongly nonlinear, cannot be obtained 
perturbatively. Rather we will find the solution to eq. (|39|) which converges to the exact one 
as e — > 0. Indeed, once e is set to zero, eq. (^9]) contains no parameters at all and should be 
solved exactly. Before we find this asymptotic {e ^ 1) solution to eq. (|39D , it is useful to 
establish some of its simple properties. 



^It should be also mentioned that after some transformations, eq.(|39D can be reduced 
to a form resembling a well known equation in radiative transfer theory, the so called 
Ambartsumyan equation which is also solved normally by iterations ( [Ambartsumyan et al. 
|1957|) . The analogy is, however, not complete and the solutions are different. 



-14- 



Consider the following function in the complex r- plane. 



G{t) = / ^ ^ (40) 

^ ' ' t' + tt'F(t') ^ ' 



Under obvious integrability conditions for the unknown function F (see e.g., [Muskhelishvili 



1953|) , the function G is a holomorphic function in the whole r- plane, cut along the 



line (— 1/e, — e) and has the degree -1 at infinity {G ~ 1/r, r — oo). This function is 
discontinuous across the cut and the jump is equal to AF = — 27ri/rF(— r). Using the 
Plemelji formulae from eg. (|40|) we thus have 

^/ N ^ f^/"^ d'T' 1 in , , 

G(-T±iO)=V — -T— — -, 41 

^ ^ Je t' -T t'F{t') tF{t) ' ^ ^ 

where, as indicated, the integral should be taken as a Cauchy principal integral. From 



eq.(|39D we infer that the function F — {Au/ Au*)^eto/ A* possesses the same analytic 
properties as G. Thus, the solution F to be found from eq. (|39|) has two branch points at 
— l/e and —e and the jump AF across {—1/e, —e). 



4.3. Solution for the spectral function. 



As we emphasized, the only small parameter available in eq.(|39D is e ^ 1. This is a very 
small parameter indeed! For shock acceleration in situations of the astrophysical interest 
the magnitude of 1/e: ~ \fpi/PQ may be as large as 10^ (for typical SNR conditions). For 
a galactic wind termination shock this parameter may be even several orders of magnitude 
larger (see [Jokipii fc Morfill 19851) . It is therefore tempting to set £ = in eq . (^9|) . In 



singularly perturbed problems, however, it is good practice to keep singular points at their 
'exact' positions in all orders of approximation. The singular character of the perturbation 
in the parameter e is quite obvious since the case e = corresponds to the absence of a 
cut-off momentum [pi = oo), and the solutions are structurally different from those for 
e 0. Such physically critical quantities as Pc and the precursor length may diverge as 
£ — i> 0. We therefore keep e small but fixed at some critical points in our further analysis. 

First, we represent eq.(p9D in the following equivalent form 

r°° dr' 1 
^'^^ ^ Is t' + t t'F{t') ^ 

+ r + ^— (42) 

Au*\ A* \Js Jl/eJ t' + tt'F{t') ^ ' 

Shifting then the lower singular point to the origin (r + e = y) we rewrite the last equation 
as follows 

Jq y' + y y'F{y') 
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where 

R 



1 



Au* ^ A* J2e F{y') y' + y-2ey' -e 

1 r°° dy' 1 
7170 y'F{y') y' + y 

F{y) = Fo{y) + ^F^ + ... (45) 



After the ansatz 
for Fq we have 

lo y' + yy'Fo{y' 
Taking ^y negative from the last equation we obtain 



f„(_,±.0)=pr^^^^ (47) 

Jo y'-yy'Fniy') yFaiy) 



y'-yy'Fo{y') yFo{y) 

Separating real and imaginary parts we get 

5^^o(-y) = V r \ (48) 

Jo y'-yy'Fo{y') 

QFoi-y±tO) = T^T^ (49) 

yFoiy) 

It is seen that the function 

(50) 

satisfies both last equations and, therefore, eq.(HB|) exactly. 




Returning to the variable r we have Fq{t) = y 7r/(r + e). This solution is not uniformly 
valid and deviates from the true solution at large r ~ e~^, i.e. the formal validity range is 
|r| <^ 1/e. The origin of this nonuniformity is obviously in the transformation from eq.(|5^ 
to eq.(0) in which only the lower singular point t = —e was kept at its correct position 
in accordance with the analytic properties, considered in the preceding subsection. The 
second singular point, namely r = —1/e has been shifted to infinity. U y ^ 1/e then the 
term ^/eR in eq.(|^) must also be taken into account. The standard way to get a uniformly 
valid expansion is to treat the solution at large r ~ 1/e in exactly the same way as we did 
at the left end of the interval e,l/e and to construct then a composite expansion. This 
requires, generally speaking, the calculation of higher order terms of the expansion, i.e. Fi 
etc., since Fq might become so small that Fq ~ y/eFi (see eq.(|^)). Our knowledge of the 
analytic properties of the solution reduces computations dramatically. Since we need just 
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a uniformly valid zero order approximation, we may omit the second term in eq.(p9D and 
obtain the following equation 

F{r) = / 51 

^ ^ Js t' + tt'F{t') ^ ' 

Here F denotes the uniformly valid zeroth order approximation to F in eq.(^). Assuming 
r <^ 1/e, from the above results we can write 



Next we observe that the transformation r i-^ 1 /r, F i— tF maps eq. (|5T|) onto itself. To be 
specific we transform (r, F) H) as follows 

e=l/r, H{i) = F/i (53) 

and for H we obtain 

rl/e (jC' 1 



According to eq.(B2|) from the last equation we infer 

^(0 = 



TT 



which is a good approximation for |^| <^ 1/e. Returning to {F,t) and taking eq.( into 
account we obtain the uniformly valid asymptotic result 




V- 55 

€){l + eT) ^ ' 

It may be seen that this method of uniformization of the asymptotic solution of eq.(P3) 



works very well in the interval(— 1/e, — e) and, what is important, restores the correct 
position of the branch point —1/e that had "escaped" to infinity in the asymptotic result 
(0). Besides that, the solution (^) exhibits a proper behavior at infinity according to 
the analytic properties discussed in the previous subsection. This method, however, might 
be still insufficient in the region r ~ 1/e, especially if t^j A" is large, and if then the ^Je 
term on the r.h.s. of eq. (|43D must be taken into account. On the other hand, the region 
r ~ +l/e corresponds to the low momenta p ~ Po and in the efficient solutions particles 
from this part of phase space play no dynamical role. We will not consider this region here. 
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5. Structure of the shock transition 



The solution for the spectral function V{p) obtained in the preceding section contains 
the undetermined constant A* that appears in the integral equation (pB]). This constant 
is in fact an eigenvalue of the nonlinear (singular) equation (|33D P|. Hence, the solution for 
V{p) found in the last section exists only if eq.(p7D possesses real solutions for A*^. Using 
the solution (pSl) and the notations (|38D we rewrite the equation for A* as follows 



(56) 



A* = A exp 




where 



i/£ dr 



tF{t) 



(57) 



This integral can be calculated conveniently by transforming it to an integral along the 
cut (eq.(^9D) and then substituting F from the approximate solution (|55D. This procedure 
yields / = \Jt^ /s. On the other hand the same integral can be calculated 'exactly' by 
multiplying eq. (|5T|) by 1/F and integrating both sides between e and 1/e. The result reads 



e 



2 



(5^ 



or / ~ yS/e. This insignificant discrepancy is de facto explained in the last paragraph of 
the preceding section. Namely, the corrections to the asymptotic result ( ^5|) in the region 
e < — r < 1/e are of the order of ^/e ^ 1. On the other hand, being integrated over an 
interval of the length 1/e, this correction produces the contribution l/\fe. This correction 
does not change the asymptotic behavior of the solution ( p5D in the dynamically most 
important part r ~ e. Substituting J = J'2./e into eq.(^) we obtain 



A" 



where 



Au* 



exp 



n 



'A* 



(59) 



(60) 



Clearly, eq.(|59|) not always has solutions and we need the shock and subshock compression 
ratios r = u\/u2 and = Uq/u2 to solve this equation for A*. 



^The fact that eq.(|33D possesses an 'inhomogeneous term' on the r.h.s. is irrelevant in 
the sense that the equation is nonlinear and that Fredholm's alternative cannot be applied. 

'^This is irrelevant for the inefficient solution, that can be directly obtained from eq. (p^) 
as a Neuman series in A provided that this parameter is sufficiently small. 
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5.1. General equation for the shock structure. The problem of uniqueness 



From eqs. (^8]) and (|55|) we have 



Jis) = \ 1 TT r 61 

^ ' V (/^os + h) {kqS + to) ^ ' 

where kqS = t. Using the relation (|3^) and inverting the integral in eq. (p8|) for F'{^) we 
obtain 



F'(^) = / e^*J(s)c/s (62) 

The last integral can be transformed to the following integral along the cut (see sec. 4. 2) 

F'(^) = — -/ e"'[J{s + tO)-J{s-iO)]ds (63) 
2m J-to/Ko 

Substituting J(s) we deduce 

F'(^) = -ae-'^+"^/o(a-^) (64) 

where 

Am* / — 1 , , 

a = JntoA*, a± = — (to ± ti) (65) 

and Jo denotes the modified Bessel function. Using the boundary condition 

/dx = Mo, = 0+, from eq.(^) we now obtain the following differential equation for the 
fiow potential ^ 

^=Uo + a e-"+* 7o(a_^')^^^' (66) 
dx Jo 

which can be always integrated in closed form. We shall describe the shock structure in the 
next subsection in more detail. Here we determine the global characteristics of the fiow. 
Using the boundary condition /dx — mi as \1/ oo and calculating the integral we 
arrive at the following relation for the compression ratios 



hr A* 

rs + (r:-l)^^ (67) 



^It is worthwhile to note that this integral diverges as if e = 0, i.e. ti = emphasizing 
the singular character of the e- perturbation. 

^Again, it is in principle possible to calculate here the integral without using the 
approximate solution (|6TD by the method applied for calculating the integral I cf. eq. (|57|) . 
The result would also differ by the factor \J2/t{. We ignore this small difference and use the 
eq.(|66D in what follows since it exactly conforms to the flow proflle. 
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where rl = + Vq/u2 — = uq/u2. It should be emphasized that we have neglected here 
H = ^(Po) in comparison with mq; not because our solution V{p) has such a property- 
it may be inaccurate at p ~ see the last paragraph in the preceding section- but 
rather by relying on an estimate that is made in Appendix C, eq.|C8|. Eg . (|67|) should be 
solved together with eq.(^) and with the subshock R-H relation (y). If the subshock is 
significantly diminished, we note that this implies automatically r ^ rg, because of our 
constraint (|^) on the Mach number v <^ 1. It should be emphasized that the condition 
r ^ Ts is a characteristic of the strong CR shock modification rather than a decrease of 
as it is often asserted in the literature. As we will see, a very strong shock modification 
may occur in the high Mach number limit also despite of ~ 4. From eqs. (|59|) and (^) we 
then obtain 



A* = r]^jTrpQpi/eexp ( -— = 1 (68) 



The numerical factor 6 here was introduced earlier in eq.( 35). It will be determined after 
the particle spectrum has been found. 

Denote il/ ^/A* = z. We then have the following bifurcation equation for z 

Tze-' = 1 (69) 

where T = 'qy/wp^p^O^ /Vt. Eg . ([69|) has either two solutions, (if F > e) or no solution at all 
(F < e, see, however, below). At the bifurcation point F = e there exists a double root. 
Therefore, these two solutions exist only if 



riV^i > B ^ (70) 

(rs - 1) ' 

where 

B = 3e/| (71) 

We implied again that Vq <^ Au which will be shown to hold unless the subshock is very 
weak. In the case F — e -C 1 this pair of solutions may be written as follows 

(72) 



l±v'2(l-e/r) 

Since r] is bounded as rg ^ 1 (MV95), a very strong subshock reduction is prohibited by 
the criterion ([70|), and, therefore, a finite subshock must remain. We shall return to the 
issue of subshock reduction later, after the fiow profile is obtained. 

This was the calculation at criticality, i.e. where F ^ e which occurs, say, at certain 
Mach number M, provided that pi is fixed. Consider now an acceleration regime that can 
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be characterized by the condition F ^ 1, i.e. far from the bifurcation point. We then have 
from eq.(^) for the first solution {z <^ 1) 

VA* ~ ri^Jirpopi/e (73) 

and 

Vl^.A (74) 



for the second one {z ^ 1). For the sake of convenience we call the solution (|73| , |72| 



"efficient" and the solution (74,72|+) "intermediate". As we argued earlier there exists also 



the third one, i.e. a conventional test particle or "inefficient" solution of eq.(|33|) at least for 
sufficiently small rj. This solution cannot be obtained by means of the approach developed 
in Sec. 4. 3 but it can be found directly from eq.(p^) as a Neuman series for A -C 1. The 
relations between these three solutions will be considered in Paper IF In the case of a 
weak influence of the CRs on the flow the inefficient solution was studied perturbatively by 
Blandford (1980) and Heavens (1983) on the basis of the diffusion-convection equation (^. 
Here we concentrate on the efficient solution for which from eqs. (pTf) and ( [73D we deduce 

r ~ Ts + 7r(rs — l)rjpi/6 (75) 

Note that the parameter po which is irrelevant for the efficient solution has indeed 
disappeared from the result, as discussed in Appendix A. The total compression ratio and, 
therefore, the acceleration efficiency depend only on the subshock compression ratio, on the 
injection rate ri{rs), and on the upper cut-off momentum pi. Combining eq.([75|) with the 
subshock R-H relation (^) we obtain an equation for the single variable w, given by 

w = — M"WT (76) 

The equation reads 

w = /i + A (l - , (77) 

where 

, = M-*«1, A= , (78) 

(7 + 1) OM— 

This equation obviously possesses only one root, that belongs to the interval (/i, 1). Since 
fi {rjpi is always large in the case of the efficient solution) we have 

w ^ /i A ^ A (79) 

in the case A -C 1. In the opposite case A ^ 1 we obtain 

w ~ (80) 

A + 1/(7 + 1) ^ ^ 
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Both expressions match at A ~ 1. Thus, the parameter A regulates the degree of shock 
modification for the efficient solution. In particular, for A > 1 one may write 



- ~ (81) 



The subshock strength can in principle be reduced strongly in this case 



3 1 

1 + 7T (82) 



' l + 3w^/^ 4 A 

We have put 7 = 5/3 in to the R-H relation (|), for short. However, A itself depends on 
through 7] but remains bounded when ^ 1. Unfortunately, the existing injection theory 
essentially implies a finite subshock, so that it is difficult to judge the behavior of 77 at small 
Ts — 1. As we shall see, rg — 1 cannot, in fact, be very small, so that the dependence ri{rs) is 
not so critical here. We may rewrite eq.(|8^) as 

r^-lc:,-——- 83 

Again one sees that formally the subshock can be quite weak for very large pi. Clearly, 
the details of acceleration in the regime rg — 1 ^ 1 may be quite sensitive to the injection 
mechanism at the weak subshock. Generally, eq.(^) places both upper and lower limits 
on the parameter pi. To see how it works consider the limit A <^ 1. In this case the total 
compression ratio is almost independent of M, whereas rg is close to four: 

r Stt 4 

-^-,p,»l, rg^^^^ (84) 

(It is implied that at least rj^/poPi > 1, see eq.(^), and therefore, rjpi ^ 1). The subshock 
strength in this case is practically not affected despite a very strong modification of the 
fiow in the precursor. If we now raise pi so that A becomes larger and rg — 1 smaller, the 
condition ([70| ) may be violated because of the large r.h.s. Denoting rj{rs = 4) = rjo, and 
rji as the value of 1] at the smallest rg to satisfy the condition ([70D, we may express this 
condition in terms of pi as follows 

, ^ 3 

1 / q\ 2 _ nz~ 

(85) 

This puts a limit on the Mach number below which the efficient-intermediate pair of 
solutions does not exist 
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It can be seen that the magnitude of injection in both the cases of modified and unmodified 
subshock is critical for the existence of efficient solutions. At the same time, the minimum 
strength of the subshock that is reached when pi approaches its upper limit in eq.(p5D, is 
not very sensitive to the parameters: 

The parameter rjipQ is presumably rather small. It characterizes the injection efficiency for 
a significantly modified subshock and preserves the subshock through quenching injection 
eventually. One sees that in practice the minimum subshock strength rsmin — 1 is of the 
order of unity since the Mach number is constrained by the condition A 1, or equivalently, 



5.2. Flow profile in the CR precursor 

The shock structure usually emerges from the balance between the nonlinear advection 
and the dissipation of the fiow energy. The dissipation is normally provided by viscosity 
or heat conduction or by both. Heat conduction alone may or may not ensure a smooth 
shock transition whereas viscosity alway does, and a viscous subshock should be inserted 
whenever other transport phenomena like heat conduction cannot provide a smooth shock 
transition. The shock amplitude and its thickness in media with quadratic nonlinearity are 
related by ul ~ z/, where z/ is a relevant dissipation coefficient. 

We started our consideration from these well known facts of shock wave theory, 
because a very similar situation arises in shocks strongly modified by CRs. For example the 
relation uil ~ ki holds in this case. In contrast to the usual gasdynamics the fiow energy is 
absorbed by the high energy particles, that have diffusivity ~ ki = n{pi). At least in the 
case of sufficiently high Mach numbers considered throughout this paper, inviscid smooth 
transitions are not possible in CR modified shocks for any finite pi. There are two reasons 
for that. First, while ^ 1 the nonlinear shock structure associated with the acceleration 
process disappears as well (r r^, eq.(^5])). Second, this structure also disappears if 
injection becomes less and less efficient as the subshock vanishes. 

First, we consider the fiow structure in a region where x is not very small Uqx/Aukq ^ 1? 
i.e. basically beyond the diffusion length of injected particles. Instead of eg . (|66|) one may 
write 

f-"«+-vf*(^) '''' 

where 

= r e'^^dt (89) 

VTT JO 
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and kI = K,{pi)/6. As in the previous subsection we focus on the efficient solution for which 
we obtain (see sec. 5.1) 



uo + {ui - no)$ [^s/-^/nlj (90) 
For ijj < kI we then obtain the following implicit solution for the flow potential \I^(a;) 

a r— , / a 



where 



a; = — V^-ln 1 + — (91) 

2uo Uq \ Uq ' 



a = 2^^ (92) 



For a\/^ /uq < 1 one obtains \1' ^ uqx + (2a/3Mo)(Moa^)^^^ although this is correct only for 
Uqx/Aukq ^ 1. The velocity proffie in this region has the following form 



u{x) = uo + 2r]Au^—^jTr— (93) 

This singular behavior may only take place in the interval 

Au^UoX^/uo_y 9 

that not necessarily exists. The prerequisite for this region to exist is a significant subshock 
reduction that is hardly possible as the preceding subsection suggests. For smaller x one 
has to use eq.(^) that yields 



u = Uq + Txri\J UqAus^ — — — (95) 

This was a thin region close to the subshock, related to the diffusion length of low-energy 
injection particles, whose dynamical role is assumed to be unimportant and no major 
change in flow speed occurs at this scale. It is perhaps worth mentioning that, if it were not 
so, the flow structure in this region would be strongly influenced by the details of injection 
which are not considered in this paper. 



For larger x, when (a;/Mo)v^' > 1, the velocity proffie is again close to linear 



a^x 



m(x) ~ a V2 + Mo 1 + In (96) 



which in the case of large Ui ^ uq and for x < I can be simplified to 



X 

u{x) ~ ui— (97) 
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where the total scale height of the precursor 

1 = ^ (98 

This velocity behavior is valid for uq/ui < x/l < 1. It is the most significant part of the 
shock transition. It comprises not only the largest velocity variation. It is also the most 
certain part regarding the asymptotic methods used here (see also below). As mentioned 
earlier, besides this longest spatial scale there also exist two short scales associated with 
low energy injection particles (cf. eq.(p3)) 



n — f^o — 2"' '2 — /^o-^ — n-. — 5 {^^ 

Far from criticality, i.e. when rj^/popi ^ 1 or when the subshock is not modified very 
strongly, the interval [h, I2) collapses and these particles play no dynamical role in the shock 
structure. Generally, this quasi- singular behavior does not seem to be persistent in the 
efficient solution. However, this conclusion might change in the case of the inefficient or the 
intermediate solution and may be useful in interpretation of observations and simulations. 
It should be recognizable as a region of very large du/dx just upstream from the subshock 
and generally should indicate that the acceleration process is close to criticality. However, 
we shall not consider this situation in more detail here. 



Beyond the precursor length, x > I from eq. (|90D we obtain 



.^.1-— — y-exp(^--j (100) 

Again, this region is less certain than that where uq/ui < x/l < 1, since for x ^ I the 
underlying asymptotic solution of the diffusion-convection equation needs some minor 
modification (see Sec. 3). At the same time m ~ Mi at this distance and the shock transition 
is virtually completed there. Furthermore, this outermost part of the precursor is accessible 
only for particles that in reality may begin to leave the system and whose behavior is 
idealized to the greatest extent by the introduction of the sharp cut-off at p = pi. We see 
that there is no major physical reason for improving formula ( p.00| ) in the present level of 
description. 



6. Particle spectrum 

The particle spectrum is the main goal in any acceleration scheme. A formal expression 
for the particle spectral slope q = SP was given already in Sec. 3 (see eq.([19|)) in terms 
of the spectral function V{p) = V{p) + Am. Its definition however involves (3 itself 
and, what is even more important, the flow structure across the whole shock transition. 
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Then an independent equation for V{p) was derived and solved, but the solution still 
contained the undetermined nonlinear eigenvalue A* that we were able to specify only 
after the shock compression ratios had been obtained. Moreover, the solution ( ]T7| ) to the 
diffusion-convection equation is strictly valid only under specific assumptions concerning 
the spectrum itself. These assumptions can thus be justified only after the spectrum is 
calculated. The close link between all these quantities is a natural aspect of nonlinear shock 
acceleration. It is this complexity of the acceleration process that has been standing in the 
way of a full analytic calculation of the particle spectrum for a long time. 

Using the same normalization for go{p) as in eq.(|26|), that is go{po) = 1, we obtain from 
eqs.(0) and (^) 

^ \ 1 [3^2 fP dp' ] 

For p ^ po we have J ~ 1, and the particle distribution is go ~ {p/po)~^^^^^^, matching 
smoothly the thermal plasma at lower momenta via the solution of the injection problem, 
as it is explained in Appendix A. The same formula ( |101|) is valid also for po < p < pi, 
while J{p) was given in Sec. 4. However, as explained, the solution for the function J{p) 
is strictly valid only when J{p) 3> 1, i.e. for sufficiently large p. We seem to even have 
some difficulty in determining the magnitude of the spectrum in the momentum region 
where J{p) 3> 1, since J must be integrated in the exponent over the region where it is not 
determined accurately, i.e. in the region where J{p) — 1 ~ 1. J can be regularly calculated 
also in this case by continuing the e— expansion in Sec. 4. At the same time this region is 
dynamically unimportant in the case of the efficient solution and there must be a way to 
get the high energy asymptotic result, i.e where J ^ 1 in eq.([101| 



without resolving this 

intermediate momentum range. In fact we already exploited this approximation solving the 
integral equation ( |55| ) for J. Namely, all necessary information about the behavior of J in 
the intermediate momentum range is 'hidden' in the eigenvalue A*, cf. eqs.(^) and ( |60D 
we may rewrite eq.( |101| ) as 



9o{p) 



J{p) 



exp 



+ 



3^2 



A* Au* 



pi 



dp' 



p'J{p') 



(102) 



and for p ^ po we may now use the solution (|6ll) for J. Far from criticality (see Sec. 5), 
when the efficient and the intermediate solutions are well separated, and focusing on the 
efficient solution, we have 

m(P) - ^ (103) 

It is seen that for the efficient solution the softening of the spectrum in the region p ^ po 
due to shock modification has no significant consequences for the high energy [p ^ po) 
spectrum; J rises so sharply that the exponential factor in eq. (|101|) remains ~ 1. Using 



- 26 - 



eqs. 



DAI 



51), and 



we then obtain from eg. (|6TD 
Jip) = 



1 + 



1 d\nj 



(104) 



3 dlnp ' pi 

We observe that for p -C pi, J ^ Cy/P satisfies this last equation; C is some constant. One 
may also solve this equation exactly, transforming it to a linear equation for J^. Indeed, 
introducing 



9' 



y 



J^, and C 



6^ 
Pi 



we obtain 



Idy 



(105) 



(106) 



After simple calculations the solution that is regular at the origin can be written down in 
the following form 

6!(-l)" 6! 



y 



6 

E 



-6C 



„=oC"(6-n)!6" 66(6 
This allows us to determine the constant 9. Since Vi ^ uq from eq.(|35| 

ld\nJ 



we have 



9 - 1 



Which can be calculated to yield 



9 = 1 + 



3 dlnp 

ICdy 
Q ydC 



p=pi 



1.09 



(107) 



(108) 



(109) 



For practical purposes, instead of eq. (|107|) , it is more convenient to use an approximate 
formula that follows immediately from eq. ( |104| ) 



J 



Pi 



P 



9 \ pi 



(110) 



where /3 ^ 1/6 for p -C pi and (3 = 9 — 1, for p ~ pi. According to eq.( |103| ) the particle 
spectrum at the shock front in this case is simply qq ~ 1/J and the full spectrum can be 
given by 



g{x,p) 



y/9 



l+(3 + 9^ 
Pi 



exp 



1 + P 
k{p) 



111) 



rrrjy'pip 

Clearly the spectrum hardens with the distance upstream and it is well localized at the 
upper momentum cut-off for ^/^(pi) ^ 1. The most interesting aspect of this solution 
is that the postshock spectrum go{p) is a power-law in a broad momentum range with 
a universal index coinciding with the test particle result for the shock of compression 
ratio 7. For p approaching pi, the spectrum Qq flattens. Nevertheless, its slope remains 
universal, decreasing from 1/2 for po ^ P ^ Pi to about 0.3 at p = pi, independently of 
any parameters. 
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7. Discussion and conclusions 

In this paper we have found that kinetic solution for diffusive shock acceleration which 
emphasizes a very strong coupling of accelerated particles with the gas flow. We have also 
demonstrated that there are three different solutions in a specified region of parameter 
space. This region is undoubtedly of great astrophysical interest. Therefore, the question, 
how efficient shock acceleration is, implies automatically the question which solution among 
these three is the attractor of a system evolving in time. This question can be also put 
in the following context. Is shock acceleration a passive process in which only a limited 
fraction of the shock energy is transformed into more or less parasitically accelerated 
particles, or is the shock itself radically restructured by these particles and they consume 
almost all its energy? Are transitions between these different states possible and if so, how 
do they occur? Thus, the notion of critical phenomena appears to be quite plausible for 
nonlinear shock acceleration. 

The first step of any quantitative theory of such phenomena should be the identification 
of a system parameter and the determination of its critical values. This parameter (or 
parameters) should encompass all the physically relevant quantities. As it was shown in 
the present paper, the most important system parameter is A = r]piM~^^^. For instance, if 
we consider the case in which M is fixed and sufficiently large, so that A < 1, acceleration 
is essentially governed by the parameter Ai = rjpi (see sec. 5). In particular, the total 
compression r oc ripi in this case. Other regions in Ai,M parameter space are discussed in 
sec. 5, and one important statement is very simple: there always exists a critical A^i\m) 
such that for Ai < A^^'' only the inefficient solution is possible. Moreover, there is an upper 
critical value for pi, or more precisely for pirjl^'^ (see eq.(^5D), beyond which the efficient 
solution fails to exist again. The reason is the following. Once this parameter exceeds the 
critical value, the backreaction of the CRs on the foreshock flow becomes so strong and, as 
a result, the spectrum at lower energies so steep, that the required CR pressure cannot be 
maintained. Also at Ai > A^^'*, the system was found to be in quite different states, ranging 
from r oc Ai, Ts ~ 4 independent of M (when A ^ 1), to r oc M^/^ with a noticeable 
subshock reduction Ts — 1 ~ 1, in the opposite extreme A > 1. 

This regime (r oc M^/^) could basically be quite natural also for the two-fluid 
model ( [Axford et al. 19T7| ; prury fc Volk 19"HI| ) since it corresponds effectively to the 
case Al = oo (in the two-fluid model it is implied that pi oo), and once the subshock 
is modified, one could very easily estimate also the total shock modification as (eq.(^) 
r oc M^^*. It is important to emphasize that such an estimate results from only two 
assumptions. The first says that the subshock is noticeably modified, and the second that it 
is not smeared out completely. No further calculations are needed. It means unambiguously 
that under these circumstances the compression ratio of a steady large Mach number shock 
must be rather large. In the two-fluid model, however, rg always crosses unity in the case of 
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efficient solution already at M ~ 10 even ii rj —* 0. [] One is forced to take that branch of 
two-fluid solution of R-H relations in which rg = 1 beyond this critical M and the scaling 
r ~ M^/^ cannot be recovered since it would correspond to a rarefaction wave rg < 1. 
In contrast, the inequality rg > 1 is always guaranteed in the kinetic description for any 
finite Ai. This is a fundamental difference between kinetic and hydrodynamic descriptions. 
The latter is not capable of preserving the kinetic link between injection and high-energy 
particles. It allows rg = 1 which means on the kinetic level of description an infinitely 
steep suprathermal particle distribution for p ^ poQ- This must normally result in Pc ^ 0. 
However, this does not occur within the two-fiuid description in which CRs and thermal 
plasma are coupled only hydrodynamically, not kinetically. In the kinetic picture presented 
above, this link is established mathematically through the nonlinear integral equation that 
simply has no eigenvalues corresponding to efficient solutions when the subshock is too 
weak and the spectrum in the suprathermal region too steep with the implication that the 
CR-pressure does not suffice. The kinetic link of suprathermal and high energy particles is 
a constitutive aspect of selfregulation in nonlinear shock acceleration. Indeed, it leads to 
the following universal (Mach number independent) relation between the flow deceleration 
in the smooth and in the discontinuous parts of the shock transition (eq. ([f5|)) 

= a^Pi (112) 

U0 — U2 tl 

Therefore, the modiflcation effect disappears whenever does the subshock. In fact this 
happens even earlier, since the criterion ( |70| ) will be violated and this solution disappears 
abruptly when uq — U2 drops below some critical value (eq.(|87D). 

One remarkable result that came about concerns the downstream spectrum go{p)- Over 
a wide portion of momentum space the power-law index depends on absolutely nothing! 
It equals 3| for the standard normalization of particle distribution. Interestingly, this 
coincides precisely with a test particle spectrum in an unmodifled shock of compression 
ratio 7. Thus, if the power-law index would be calculated using the ratio of speciflc heats 
7c (which is clearly very close to 4/3) rather than the compression ratio, one would obtain 
essentially the same result. Our solutions imply that r can be large, at least larger than 7. 
Nevertheless the spectrum fails to flatten with growing compression. The resolution of this 
'paradox' lies obviously in the thickness of the shock, that grows simultaneously with r and 
Pi, and these two oppositely acting factors compensate each other exactly. Similarly there 



''It is noteworthy that the sequence of limiting transitions is important here, i.e. in the 
process pi — 00, — the parameter Ai remains infinite. We recall in this regard that the 
problem under question is a singularly perturbed one in the parameter 1/pi. 

^This argument is not correct when the gas fiow develops a very large gradient just in 
front of the subshock due to the pressure of suprathermal particles. Similar situation was 
briefly considered in Sec. 5. 2. 



-29- 



is no flattening of the spectrum with momentum unless p ~ Pi, since particles that sample 
larger compression at larger momenta see a more extended flow structure as well. One could 
say that while in the test particle theory the energy spectrum is scale invariant because 
all the internal shock scales are irrelevant, in the nonlinear theory both the momentum 
spectrum and the shock structure are scale invariant in a certain region of phase space. The 
spectrum flattening at the cut-off {q decreases from 3.5 to about 3.3) is obviously caused 
by the presence of this sharp cut-off itself. The mathematical nature of this flattening is 
very simple and can be explained using the analytic properties of the spectral function 
J{p) and hence, those of go{p) in the complex p-plane. Namely, in an unphysical region 
at p = —Pi, do has the branch point, which is formally responsible for the flattening at 
p — Pi- In a physically more realistic case in which the losses are distributed at p ~ pi, this 
flattening may not appear. (We therefore do not speculate for example about a possible 
link of this spectrum flattening with the bump around the energy lO^^eK on the galactic 
CR spectrum.) In the present picture, however, particles with p ^ pi seem to 'know' how 
high is the cut-off, which appears to be strange from the point of view of causality. This 
'information' is provided again, by the flow structure, in fact through the singularity of 
the spectral function V{p) at p = — pi. As particles come nearer to pi they also diffuse 
to the outer region of the precursor where u{x) approaches Ui exponentially. Hence these 
particles develop a spectrum that is closer to g = 3r/(r — 1). Recall that in the case under 
consideration r ^ 1. Moreover, if we formally extend go{p) beyond p = pi, considering go as 
a test particle spectrum, we get the power-law index g = 3 as p — > oo instead of the 'exact' 
g = 3/(1 — r~^) ~ 3. This small deviation can be easily removed and the modiflcation 
needed for that will not affect the region where q = 3^/2- This difference in the spectrum 
slope occurred because, for simplicity, we treated the dynamically unimportant region 
^ > Ki on the same basis as the region ^ < ki (see the end of sec. 5). 

It would be interesting to examine what may happen in a time dependent case. 
Formally, the generalization of our results to time dependent acceleration should not 
be a serious problem. Namely, for sufficiently large pi one can adopt an 'adiabatic' 
approximation in which the sohition has approximately the same form as that found in this 
paper, except for parameters this sohition depends on; they will slowly evolve with time. 
If the Mach number can be fixed, the only parameter we need is pi(t), all the others and 
the flow conflguration are calculable. The cut-off increment may, at least tentatively, be 
obtained from the usual equation dpi/dt — Pi/Tacc(pi)- Since the acceleration time equals 
Tacc ^ i^{pi)/u^i the process should be slow compared with the time scale at lower momenta. 
This should justify an adiabatic approach. However, the spectrum at the cut-off may differ 
from our steady state solution significantly. 

Our main conclusion concerning the overall spectrum to have a universal index q = 3^/2 
may also be a subject for corrections due to a number of physical reasons ignored in our 
consideration. These should be taken into account before we attempt to compare our results 
with any observations or numerical calculations. We have studied the simplest case of a 
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plane steady shock with specified Bohm diffusion in which the magnetic field either plays 
no dynamical role or is quasiparallel to the shock normal. Clearly a more realistic geometry, 
time dependence and losses increasing continuously with the particle energy, will steepen 
the spectrum. We also have simplified our description of subrelativistic particles replacing 
p{l by unity in the pressure integral (eq.|26|) , (see eq.(p3D and the text below it). 

Doing so we could overlook some feature in the particle spectrum at p ~ 1 and ignore the 
role of a possible relativistic peak at p ^ 1 in the partial pressure. Generally, this peak may 
or may not develop in a modified shock. In the cases considered throughout this paper its 
contribution to the particle pressure would be much smaller than that of the high energy 



particles p <^ Pi (see Appendix A in [Malkov fc Volk 1996|) . The most probable correction 



to the particle spectrum due to the relativistic peak should be a turn in the spectral slope 
at p ~ 1 just because of the correspondent factors in the kernel of the integral equation 



(0). And finally, particle diffusion is normally assumed to be provided by waves that, in 
turn, are excited by accelerated particles. This link may not only soften the spectrum but 
also change other characteristics of acceleration process noticeably. Notwithstanding the 
physical importance of the above factors, we believe that the key issue for understanding 
the nonlinear shock acceleration is the gas flow-CR coupling. This coupling may, at least 
formally, operate under prescribed injection, diffusivity and upper cut-off momentum. It 
should be understood as such before a comprehensive nonlinear theory of shock acceleration 
is done. The major reason for that is of course the multiplicity of the solution. 

We have presented only the steady state solutions and specified the conditions under 
which they appear multiply. The main issue now is of course their realisability. Whereas 
for the case of a unique solution a proof of its stability suffices, in a parameter range where 
all the three solutions are admitted, the full time dependent calculations may be needed to 
predict the system behavior. 

Already at this stage several important features of the acceleration process can be 
pointed out. First, there are no such 'minor' things like the injection rate or the cut-off 
momentum that can be considered as having a small effect on the acceleration process 
and their treatment should not be oversimplified just by arguing that the former is small 
and the latter is large. These parameters enter all the main results and in many cases 
symmetrically, e.g., through Ai = rjpi. Therefore, both rj and pi are equally important. 
This conclusion is strictly valid for a steady state. One can imagine a time evolution in 
which e.g., pi increases but t] decreases leaving a parameter such as Ai that governs the 
steady state fiow structure (or its equivalent in a time dependent formulation, which would 
have the meaning of the energy content of accelerated particles) approximately constant. 
Furthermore, in time dependent situations the constraint (|112| ) is strictly speaking not 
valid, and a complete smoothing of the subshock cannot be excluded. This last remark 
might be important for the interpretation of numerical results. 

A number of time dependent numerical solutions were obtained in the past {e.g., |Bell| 
nSTt [b'alle fc Giddings 19871 ). Except perhaps the case considered by Bell, (M = 100) it 
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is not clear whether they are within the parameter range where multiple solutions appear 
or should rather be regarded as inefficient solutions, sometimes with quite a strong shock 
modification and high acceleration efficiency. Furthermore, at moderate Mach numbers 
and near the bifurcation point the 'inefficient' solution may be not very different from two 
solutions that we termed efficient and intermediate. We will discuss the relations between 
these three solutions in more detail in Paper II. 

I wish to express special thanks to Heinz Volk for drawing my attention to the 
problem of nonlinear shock acceleration, for many fruitful discussions of this subject, for 
his thoughtful reading of this manuscript and for his invaluable comments and suggestions 
concerning this paper. This work was done within the Sonderforschungsbereich 328, 
"Entwicklung von Galaxien" of the Deutsche Forschungsgemeinschaft (DFG). 



A. Injection 

The injection rate calculated in (MV95, M96, MV97), essentially depends on details of 
the subshock dissipation mechanism like the level and the spectrum of underlying plasma 
turbulence. The obliqueness of the magnetic field and the presence of reflected particles 
( pcholer 1990| ) will clearly influence the injection rate too. These parameters may vary 
from model to model and it is therefore perfectly reasonable to consider the injection rate 
rj as another free parameter along with M and pi, at least in the context of nonlinear 
shock acceleration studied here. The subshock compression ratio rs{r], M,pi) obtained 
in this paper may then be combined with the injection rate rj = ri{rs) obtained in the 
papers listed above to yield both = rs{M,pi) and r] = ri{M,pi). However, a microscopic 
subshock dissipation mechanism that operates quasi-independently of the large scale shock 
modiflcation, should be identifled for this purpose in each particular physical model (see 
M96 for an example of such a mechanism). 

Besides the injection rate t] we introduced the injection momentum pq (eq.(^) in order 
to separate the CR and the thermal populations. This can be done unambiguously if the 
injection is weak so that the injected particles which are between these two populations 
in the momentum space, do not modify signiflcantly the flow upstream by themselves 
but only after being accelerated to very high energies and therefore over a much larger 
scale. This is the essence of our approach to distinguish between thermal and high-energy 
populations based on the difference in diffusion length of the respective components. 
Clearly the monotonic growth of the diffusion coefficient with the particle momentum is 
very critical here. If the subshock is sufficiently strong, these two populations of particles 
do not even overlap in the upstream phase space. Fig. 2. They are separated by the region 
■^^0 < Ip/^^I < a(Mo — U2), where a 1 (M96), p is measured in the local upstream 
frame where u{x) ~ Uq and Vq is the thermal velocity at the same point. Therefore, the 
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momentum pq is a rather formal lower boundary separating the CR gas from thermal 
plasma and from the viewpoint of the upstream distribution, pq (or may be Poi^), where 
6 is the pitch-angle) can be chosen arbitrarily within this empty part of the phase space. 
Turning to the downstream region we note that the injection solution yields the particle 
spectrum that starts from the thermal Maxwellian and evolves at higher momenta into 
the power-law determined by the subshock compression. This allows the smooth matching 
of this injection solution with the low- momenta asymptotic solution of eq.(^ within an 
extended overlapping interval in the region p ^ am{uo — U2), (see also M96, MV97). This 
means that again, one should be able to choose Pq arbitrarily within the overlapping interval 
and all physically meaningful results must be independent of Pq. If po indeed hes in this 
empty phase space region upstream, then Pjnj = in eqs.([5|,|26|). However, such a choice 
of Po may cause technical difficulties since the diffusion-convection equation is not valid 
in this region due to anisotropy of pitch- angle distribution. If we take po larger, just to 
make the diffusion-convection equation valid, then Pinj 7^ and must generally speaking be 
retained in Bernoulli's integral, eqs.(|^,^). It is also worth noting that the matching of the 
downstream injection solution onto the standard power-law occurs typically at p somewhat 
larger than the void in the upstream distribution. Thus, po may be slightly higher than the 
lower boundary of the upstream CR distribution. 

The above details of the phase space geometry, being essential for injection and for 
its link with the further acceleration, are not so critical for the subject of this paper, i.e. 
for the shock modification by high-energy particles. It should be emphasized, however, 
that an important physical condition for legitimating both the injection theory and the 
present study is that the suprathermal particles with momenta p ~ po, play no considerable 
dynamical role in the precursor, i.e. their contribution to Pc is much smaller than that of the 
rest of momentum space Po ^ p < pi- As we have stressed already, they do not produce a 
noticeable fiow variation over their own diffusion length k{pq)/uq which for example justifies 



^This statement should not be taken literally. First of all, po formally enters the definition 
of the CR number density ric (sec. 4) and thus the parameter 77. But this is an implicit 
dependence on pq through Uc and once ric is defined po -invariant, {e.g., by subtraction of 
the thermal Maxwellian from the total particle distribution), this dependence disappears. 
Second, the "injection momentum" Pq can appear in some criticality conditions (Sec. 5) and 
also in situations when the approximation described just below is only marginally correct. 
This means that the description of the injection process by only one parameter rj is insufficient 
in such a situations. On the other hand in selfconsistent injection models both the parameters 
rj and Pq should be explicitly specified in terms of physical quantities like the thermal velocity 
and the Mach number. If so, a complete solution of the acceleration problem depends not on 
parameters like 1] or pq but only on these physical quantities. At the same time, as it may be 
seen from Fig. 2, the transition region between the thermal Maxwellian and the high-energy 
power-law is relatively narrow, so that the quantity po has also a clear physical meaning. 
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the adoption of the injection calculations performed in a quasi-homogeneous upstream flow. 

The injection theory produces a suprathermal downstream asymptotics go{p) — Qp~'^° 
where go = 3m2/(mo ~ "^2) and Q depends also on go along with a number of other plasma 
parameters near the subshock like the downstream temperature and the magnitude of 
MHD-turbulence (see MV95, M96). This is essentially the same as eg . (|2^) . The theory of 
injection also determines the spatial distribution of g{x, p) at the subshock for particles 
with momenta k{p)/uq ~ /inj = i^{po)/uo <^ /, where u{l) ^ Ui. Therefore /mj is virtually 
irrelevant for our further consideration and only the parameter Q is needed. 



B. Diffusion-convection equation 

We examine the idea that the result (0) that is strictly valid for small values of "^/k, 
is also a good approximation to the solution of the convection diffusion equation in general. 
To this purpose, it is natural to seek the solution of eq.@ in the following form 

^ = ^o(p)exp|-^^^|^(p,^) (Bl) 

Clearly we have g{p, 0) = 1 as a boundary condition. Our objective here is to demonstrate 
that g{p, \E') ~ 1 also for \E'/k noticeably larger than unity. A practically important region 
is of course that where \E'/fi; is not very large due to the exponential factor in eq. ([Bl|) . 
Introducing also 

for the function g we have the following equation that can be easily derived from eq.@ 

where r = 3 Inp/po. At this point we need some more information about the function Q{'^). 
As we have shown in sec. 5. 2., for not very small more precisely for 

^//co > {uo/Auf/7]^popi, (B3) 

and in the case of strong shock modification, the velocity of the flow behaves as follows (see 
eq-O) 

where $ is the probability integral. Thus, for \E' ^ i^{pi), to a good approximation we may 
write 

QW ^ ^ (B5) 
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For larger \1/ > ki the function falls off exponentially. But this region is dynamically 

unimportant {g is exponentially small already for p <^ pi) and we confine our consideration 



to the case in which eg . ( p5|) holds. Introducing a new variable z = y2\l//fi: and performing 
the variable transformation (^I^,t) {z,t) eg. (P2|) rewrites 



^ + a{T)z^ + 5{T)z^g--^=0, ^(r,0) = l (B6) 



Or dz dz"^ 

where 



a=^ + 2/3(r) and 5 = (1 + - /?) - 1^ (B7) 
6 b 2 or 



Physically, the problem (|B6|) should be posed as follows. At r ~ ~ po), g{.z, r) has to 
be matched smoothly onto the solution of the injection problem as it was discussed earlier. 
Then g and hence g can be determined from eg. (p6|) provided that (3{t) is known. This 
latter, in turn, depends on the spectral index of the downstream distribution to be found 
from eg.(|12]). As it is shown in Appendix C, unless p ~ po or P ~ Pi the function (3{t) and 
thus 5{t) is a slow function of r. We observe that a characteristic 'time scale' in eg. (P^) is 
r ~ 1. Thus, taking a and 5 as approximately constant one can easily solve the problem 
( p6D in terms of parabolic cylinder functions. But even this simple solution is not really 



needed for our purposes. It turns out (Appendix C) that 5 ^ 1 for po ^ P ^ Pi- Consider 
the region r > 1, where the solution 'forgets' the 'initial' condition g{T ~ 0, 2), that is given 
by the high energy end of the injection solution. The solution to eg. (|B6|) will be determined 
by the boundary condition ^(r, 0) = 1 and one may see that for not too large z < 1/^/5 we 
also have g{T, z) ~ 1. Furthermore, the injection solution yields the spatial distribution of 
the particle spectrum that is consistent with the exponential decay oc exp(— (1 + 
(see MV95). We may thus put ^(0,2;) = 1 which yields g{T,z) ^ 1 for small 6. As we 
mentioned above this solution is not valid for z 1/ \/6, but once 6 is small we can ignore 
this region since g{p,z) ~ (?o exp(— 1/25) there. 



Finally, we note that this simple solution (^ = 1 ) of eg. ( [B6| ) is exact if we reguire 5 = 
in eg. (p7|) . This yields for the spectral index 

^ o 1 + a + (6ft - Ije-'"-" 
« = ^"W=\(l + ,3.)-(6/3„-l)e-/3 

where 3/5o is the spectral index at r ~ 0, {p — po), which is given by the injection theory and 
coincides with the standard power-law index calculated for the subshock compression ratio. 
For larger r, jS drops to 1/6 rapidly. Clearly, for ^ = 1 to hold, the above behavior of (3{t) 
must be consistent with the solution of eg.(|l2[), which will be confirmed in Appendix C. 



C. Verifying assumptions 



It has been understood for a long time that in contrast to the test particle theory, 
the shock acceleration in the nonlinear regime cannot be decomposed as a seguence of 
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independent processes like the formation of the overall flow structure, particle spectrum, 
and injection, or losses. Even the magnitude of the cut-off momentum cannot be specified 
prior to obtaining the full solution when it is influenced by the solution itself (we did not 
consider this case here). All these factors are linked very tightly in the solution obtained in 
the preceding sections and non of them can be prescribed a priori. Nevertheless, we made 
several assumptions concerning mostly the flow profile and the particle spectrum. For the 
sake of convenience we list here these assumptions along with the confirming results. 

First, we neglect the adiabatic gas compression term in the Bernoulli's equation 
upstream, eq.@. Denoting 



uix] 



Ui 



and J^{x) = 1 



{Pc + Pi 



(CI) 



we rewrite eq.(p[) as 



F{v) = v + 



1 

1 



(C2) 



We must specify the condition under which we may invert F{v) to obtain t> ~ JF. First of 



all we observe that the l.h.s. of eg. (|C2|) has a very sharp minimum aX v = Vm = M ^/('>'+^). 



Since f m -C 1 the approximation f ~ JF may be good over most of the precursor, where 
w > f m (see Fig. 3). Indeed, using eq.(p5|) we can rewrite this inequality as 



Txrj 



PoPi 



-UqX 



Kq 



> M ^+1 



(C3) 



In the trivial case (mi/mo)M"2/(^+^) < 1, i.e. uui/uq = 2M"2(ui/uo)^+^ < 1 the 
approximation f ^ JF is correct for all x > 0. If only the weaker condition z/ <^ 1 is fulfilled, 
we obtain the following constraint (eqs. (|75|JC3|) ) 



UqX 
Kq 



> 



(C4) 



Substituting pimax and rsmin from Sec. 5.1, we obtain 



UqX 
Kq 



1 

' 2 r-^ 



i)v^ 



(C5) 



That means that we may neglect the adiabatic compression term already inside the scale 
height of injected particles ~ Kq/uq. Smaller x belong virtually to the internal subshock 
structure which cannot be resolved within the diffusion-convection description in any case. 
If the nonrelativistic correction to the diffusion coefficient is important, one should write 
^0 = 1^0 /Po, where kq is the 'true' diffusion coefficient at p = Pq. Then the inequahty (|C5|) 
becomes more stringent and x must be outside the diffusion length of low energy injection 
particles. 
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Our next assumption concerned the functions (3{t) and S{t) that enter eq.( 
According to Sec. 6 and for sufficiently large p ^ we may write P ~ {p/3){d\nJ/dp), to 
yield 

^ ~ 6 ' 1 + 60^/7^1 ^^^^ 

One sees that P is close to | for p not too close to the upper cut-off. What is even more 
important for our treatment of eg. (p6|) in Appendix B is that 5^1. Even in the worst 
case possible, namely when p pi, we still have 6 ^ 0.1, so that the result ([T7|) , based 
largely on the above behavior of parameters f3 and 6, is still a good approximation to the 
solution of the diffusion-convection equation. 

A number of our approximations depend on the value Vq = V{po), Sec. 3.1. 
Unfortunately it is impossible to determine this quantity directly from our asymptotic 
solution for V{p) since it is formally valid for V{p) 3> Vq. We therefore should estimate Vq 
directly from the definition (pO]). 

^ (C7) 
=0 "0 



du 

\/o ~ — 

ax 



The quantity kq = PqKq accounts of the nonrelativistic correction to the diffusion coefficient 



used. According to eq. (|C3|) we thus have 



Vq I 

7r?7^/p^]9oVl - V'^s (C8) 

Uq " 

Since the criticality parameter Tj^/p^pl > 1 cannot be very large and Au not very small (see 
the end of Sec.5.1), we infer that Vq < Au is a reasonable approximation for sufficiently 
small pq. 
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Fig. 1. — Biffurcation curve for three possible solutions in the {R,ri) plane. Here R is the 
flow compression upstream, between the infinity and the subshock, R = Ui/uq and t] is the 
injection rate. 

Fig. 2. — Omnidirectional injection spectrum downstream (measured in the shock frame) 
calculated for the downstream temperature T2 = 2 ■ lO^K and compression ratio rg ~ 4 
(MV97). Dashed lines show the correspondent upstream distribution schematically. A 
possible region for matching is indicated by the box. 



Fig. 3. — The l.h.s of eq.( |C2D as a function of v for M = 10. 
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